#ifndef PCPS_FULLCI_HEADER
#define PCPS_FULLCI_HEADER

#include <src/pcps.h>
#include <src/matrix_actor.h>

namespace pcps { namespace fci {

  //-----------------------------------------------------------------------------------
  // pcps::fci::fullci -- Computes the minimum Hamiltonian eigenvalue and eigenvector
  //                      in the full configuration interaction basis.
  //-----------------------------------------------------------------------------------

  template <class SCALAR> void fullci(const pcps::Input & userinp, const pcps::Hamiltonian<SCALAR> & ham);

  //-----------------------------------------------------------------------------------
  // pcps::fci::HamiltonianActor -- Class for acting the hamiltonian on a trial vector.
  //-----------------------------------------------------------------------------------

  template <class SCALAR> class HamiltonianActor : public pcps::MatrixActor<SCALAR> {

    private:

      // data members
      const int _n_a_sites;
      const int _n_b_sites;
      const int _n_a_config;
      const int _n_b_config;
      const int _n_a_particles;
      const int _n_b_particles;
      const double _thresh;
      std::vector<int> _config_a;
      std::vector<int> _config_b;
      std::vector<int> _index_lookup_a;
      std::vector<int> _index_lookup_b;
      std::vector<int> _rl_ops_a;
      std::vector<int> _rl_ops_b;
      std::vector<int> _rl_ord_a;
      std::vector<int> _rl_ord_b;
      std::vector<int> _nn_ops_ab;
      std::vector<SCALAR> _rl_coeff_a;
      std::vector<SCALAR> _rl_coeff_b;
      std::vector<SCALAR> _nn_coeff_ab;

      // function to get the compound index of a given configuration
      int get_index(const std::vector<int> & ca, const std::vector<int> & cb) const {

        // get compound index of alpha configuration
        int ia = 0;
        for (int i = 0; i < _n_a_sites; i++)
          ia = 2*ia + ca[i];
        assert( _index_lookup_a[ia] >= 0 );
        assert( _index_lookup_a[ia] < _n_a_config );

        // get compound index of beta configuration
        int ib = 0;
        for (int i = 0; i < _n_b_sites; i++)
          ib = 2*ib + cb[i];
        assert( _index_lookup_b[ib] >= 0 );
        assert( _index_lookup_b[ib] < _n_b_config );

        // return combined compound index
        return _n_a_config * _index_lookup_b[ib] + _index_lookup_a[ia];

      }

    public:

      // constructor
      HamiltonianActor(const pcps::Input & userinp, const pcps::Hamiltonian<SCALAR> & ham);

      // function to return the total number of configurations
      int total_configs() const { return _n_a_config * _n_b_config; }

      // function to perform y = A d
      void operate_by_A(const int n, SCALAR * const d, SCALAR * const y, const int myrank, const int tid);

      // function to perform y = M d
      void operate_by_M(const int n, SCALAR * const d, SCALAR * const y, const int myrank, const int tid);

      // function to perform y = M^(-1) d
      void operate_by_M_inv(const int n, SCALAR * const d, SCALAR * const y, const int myrank, const int tid);

      // function to check if the iterations have converged
      // result should be computed using the root process/thread and returned on all processes/threads
      bool converged(const double residual, const int myrank, const int tid);

      // function to print a supplied vector
      void print_vector(const int n, const SCALAR * const x) const;
  
  };

} // end namespace fci

} // end namespace pcps

#endif
